home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / gpplib22.zoo / libsrc / xcomplex.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-11-30  |  6.0 KB  |  257 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17.  
  18. #ifdef __GNUG__
  19. #pragma implementation
  20. #endif
  21. #include <xcomplex.h>
  22. #include <std.h>
  23. #include <builtin.h>
  24.  
  25. // error handling
  26.  
  27. void default_Complex_error_handler(const char* msg)
  28. {
  29.   cerr << "Fatal Complex arithmetic error. " << msg << "\n";
  30.   exit(1);
  31. }
  32.  
  33. one_arg_error_handler_t Complex_error_handler = default_Complex_error_handler;
  34.  
  35. one_arg_error_handler_t set_Complex_error_handler(one_arg_error_handler_t f)
  36. {
  37.   one_arg_error_handler_t old = Complex_error_handler;
  38.   Complex_error_handler = f;
  39.   return old;
  40. }
  41.  
  42. void  Complex::error(const char* msg) const
  43. {
  44.   (*Complex_error_handler)(msg);
  45. }
  46.  
  47. /* from romine@xagsun.epm.ornl.gov */
  48. Complex /* const */ operator / (const Complex& x, const Complex& y)
  49. {
  50.   double den = fabs(y.real()) + fabs(y.imag());
  51.   if (den == 0.0) x.error ("Attempted division by zero.");
  52.   double xrden = x.real() / den;
  53.   double xiden = x.imag() / den;
  54.   double yrden = y.real() / den;
  55.   double yiden = y.imag() / den;
  56.   double nrm   = yrden * yrden + yiden * yiden;
  57.   return Complex((xrden * yrden + xiden * yiden) / nrm,
  58.                  (xiden * yrden - xrden * yiden) / nrm);
  59. }
  60.  
  61. Complex& Complex::operator /= (const Complex& y)
  62. {
  63.   double den = fabs(y.real()) + fabs(y.imag());
  64.   if (den == 0.0) error ("Attempted division by zero.");
  65.   double xrden = re / den;
  66.   double xiden = im / den;
  67.   double yrden = y.real() / den;
  68.   double yiden = y.imag() / den;
  69.   double nrm   = yrden * yrden + yiden * yiden;
  70.   re = (xrden * yrden + xiden * yiden) / nrm;
  71.   im = (xiden * yrden - xrden * yiden) / nrm;
  72.   return *this;
  73. }
  74.  
  75. Complex /* const */ operator / (double x, const Complex& y)
  76. {
  77.   double den = norm(y);
  78.   if (den == 0.0) y.error ("Attempted division by zero.");
  79.   return Complex((x * y.real()) / den, -(x * y.imag()) / den);
  80. }
  81.  
  82. Complex /* const */ operator / (const Complex& x, double y)
  83. {
  84.   if (y == 0.0) x.error ("Attempted division by zero.");
  85.   return Complex(x.real() / y, x.imag() / y);
  86. }
  87.  
  88.  
  89. Complex& Complex::operator /= (double y)
  90. {
  91.   if (y == 0.0) error ("Attempted division by zero.");
  92.   re /= y;  im /= y;
  93.   return *this;
  94. }
  95.  
  96.  
  97. Complex /* const */ exp(const Complex& x)
  98. {
  99.   double r = exp(x.real());
  100.   return Complex(r * cos(x.imag()), 
  101.                  r * sin(x.imag()));
  102. }
  103.  
  104. Complex /* const */ cosh(const Complex& x)
  105. {
  106.   return Complex(cos(x.imag()) * cosh(x.real()), 
  107.                  sin(x.imag()) * sinh(x.real()));
  108. }
  109.  
  110. Complex /* const */ sinh(const Complex& x)
  111. {
  112.   return Complex(cos(x.imag()) * sinh(x.real()), 
  113.                  sin(x.imag()) * cosh(x.real()));
  114. }
  115.  
  116. Complex /* const */ cos(const Complex& x)
  117. {
  118.   return Complex(cos(x.real()) * cosh(x.imag()), 
  119.                  -sin(x.real()) * sinh(x.imag()));
  120. }
  121.  
  122. Complex /* const */ sin(const Complex& x)
  123. {
  124.   return Complex(sin(x.real()) * cosh(x.imag()), 
  125.                  cos(x.real()) * sinh(x.imag()));
  126. }
  127.  
  128. Complex /* const */ log(const Complex& x)
  129. {
  130.   double h = hypot(x.real(), x.imag());
  131.   if (h <= 0.0) x.error("attempted log of zero magnitude number.");
  132.   return Complex(log(h), atan2(x.imag(), x.real()));
  133. }
  134.  
  135. // Corrections based on reports from: thc@cs.brown.edu & saito@sdr.slb.com
  136. Complex /* const */ pow(const Complex& x, const Complex& p)
  137. {
  138.   double h = hypot(x.real(), x.imag());
  139.   if (h <= 0.0) x.error("attempted power of zero magnitude number.");
  140.  
  141.   double a = atan2(x.imag(), x.real());
  142.   double lr = pow(h, p.real());
  143.   double li = p.real() * a;
  144.   if (p.imag() != 0.0)
  145.   {
  146.     lr /= exp(p.imag() * a);
  147.     li += p.imag() * log(h);
  148.   }
  149.   return Complex(lr * cos(li), lr * sin(li));
  150. }
  151.  
  152. Complex /* const */ pow(const Complex& x, double p)
  153. {
  154.   double h = hypot(x.real(), x.imag());
  155.   if (h <= 0.0) x.error("attempted power of zero magnitude number.");
  156.   double lr = pow(h, p);
  157.   double a = atan2(x.imag(), x.real());
  158.   double li = p * a;
  159.   return Complex(lr * cos(li), lr * sin(li));
  160. }
  161.  
  162.  
  163. Complex /* const */ sqrt(const Complex& x)
  164. {
  165.   if (x.real() == 0.0 && x.imag() == 0.0)
  166.     return Complex(0.0, 0.0);
  167.   else
  168.   {
  169.     double s = sqrt((fabs(x.real()) + hypot(x.real(), x.imag())) * 0.5);
  170.     double d = (x.imag() / s) * 0.5;
  171.     if (x.real() > 0.0)
  172.       return Complex(s, d);
  173.     else if (x.imag() >= 0.0)
  174.       return Complex(d, s);
  175.     else
  176.       return Complex(-d, -s);
  177.   }
  178. }
  179.  
  180.  
  181. Complex /* const */ pow(const Complex& x, int p)
  182. {
  183.   if (p == 0)
  184.     return Complex(1.0, 0.0);
  185.   else if (x == 0.0)
  186.     return Complex(0.0, 0.0);
  187.   else
  188.   {
  189.     Complex res(1.0, 0.0);
  190.     Complex b = x;
  191.     if (p < 0)
  192.     {
  193.       p = -p;
  194.       b = 1.0 / b;
  195.     }
  196.     for(;;)
  197.     {
  198.       if (p & 1)
  199.         res *= b;
  200.       if ((p >>= 1) == 0)
  201.         return res;
  202.       else
  203.         b *= b;
  204.     }
  205.   }
  206. }
  207.  
  208. ostream& operator << (ostream& s, const Complex& x)
  209. {
  210.   return s << "(" << x.real() << ", " << x.imag() << ")" ;
  211. }
  212.  
  213. istream& operator >> (istream& s, Complex& x)
  214. {
  215. #ifdef _OLD_STREAMS
  216.   if (!s.good())
  217.   {
  218.     return s;
  219.   }
  220. #else
  221.   if (!s.ipfx(0))
  222.   {
  223.     s.clear(ios::failbit|s.rdstate()); // Redundant if using GNU iostreams.
  224.     return s;
  225.   }
  226. #endif
  227.   double r, i;
  228.   char ch;
  229.   s >> ws;
  230.   s.get(ch);
  231.   if (ch == '(')
  232.   {
  233.     s >> r;
  234.     s >> ws;
  235.     s.get(ch);
  236.     if (ch == ',')
  237.     {
  238.       s >> i;
  239.       s >> ws;
  240.       s .get(ch);
  241.     }
  242.     else
  243.       i = 0;
  244.     if (ch != ')')
  245.       s.clear(ios::failbit);
  246.   }
  247.   else
  248.   {
  249.     s.putback(ch);
  250.     s >> r;
  251.     i = 0;
  252.   }
  253.   x = Complex(r, i);
  254.   return s;
  255. }
  256.  
  257.